home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / listings / v_10_05 / 1005044a < prev    next >
Encoding:
Text File  |  1992-03-07  |  3.8 KB  |  148 lines
  1.                    /* Listing 4 */
  2.  
  3. /*****************************************************
  4.            File Name: GS_JRD_C.C
  5.        Expanded Name: gauss-jordan fortran style
  6. Global Function List: gauss_jordan
  7.          Portability: Standard C
  8. ******************************************************/
  9.  
  10. /* Standard Library */
  11. #include <float.h>
  12. #include <math.h>
  13. #include <stdlib.h>
  14.  
  15. /* Types and Prototypes */
  16. #include <matrix_t.h>
  17.  
  18. /* Own */
  19. #include <gs_jrd.h>
  20.  
  21. /*****************************************************
  22.        Name: gauss_jordan
  23.  Parameters: A matrix of coeficients
  24.              n number of rows, cols in A
  25.              B matrix of constant vectors
  26.              m number of vectors, cols, in B
  27.      Return: SUCCESS, 0, or FAIL, -1, directly
  28.              Inverse of A in A indirectly
  29.              Solution vectors X in B indirectly
  30. Description: Gauss-Jordan elimination with partial
  31.              pivoting is used to invert a matrix and
  32.              solve a system of linear equations.
  33.              The equations are expressed in the matrix
  34.              form as [A][X] = [B].  The A matrix is an
  35.              n by n matrix.  The B matrix is an n by m
  36.              matrix.  B is replaced with the solution
  37.              vectors X.  A is replaced with the
  38.              inverse of A.
  39. *****************************************************/
  40. int gauss_jordan( matrix_t A, size_t n,
  41.       matrix_t B, size_t m )
  42.    {
  43.  
  44.    size_t i, j, k, PivotRow;
  45.    double Pivot, *Ai, *Aj;
  46.  
  47.    /* Loop through all the rows. */
  48.    for ( i = 0; i < n; i++ )
  49.       {
  50.  
  51.       /* Set the default pivot to the
  52.       ** current row. */
  53.       PivotRow = i;
  54.       Pivot = A[i][i];
  55.  
  56.       /* Loop through the rows below the current
  57.       ** row Looking for a large diagonal to
  58.       ** pivot */
  59.       for ( j = i + 1; j < n; j++ )
  60.          {
  61.          /* Find the pivot row */
  62.          if ( fabs( A[j][i] ) >= fabs( Pivot ) )
  63.             {
  64.             PivotRow = j;
  65.             Pivot = A[j][i];
  66.             }   /* if  */
  67.          }   /* for j */
  68.  
  69.       /* Check to make sure we are not pivoting
  70.       ** around the current row */
  71.       if ( PivotRow != i )
  72.          {
  73.          /* Swap the current row with the
  74.          ** pivot row */
  75.  
  76.          /* Swap the pointers to the rows. */
  77.          Ai = A[i];
  78.          A[i] = A[PivotRow];
  79.          A[PivotRow] = Ai;
  80.  
  81.  
  82.          /* Swap the pointer to the rows. */
  83.          Ai = B[i];
  84.          B[i] = B[PivotRow];
  85.          B[PivotRow] = Ai;
  86.  
  87.          }   /* if PivotRow */
  88.  
  89.       /* Check for a divide by zero error */
  90.       if ( fabs( Pivot ) <= DBL_MIN )
  91.          {
  92.          return ( -1 );
  93.          }   /* if  */
  94.  
  95.       /* Invert the pivot so we can use * operator
  96.       ** instead of / operator */
  97.       Pivot = 1.0 / Pivot;
  98.  
  99.       /* Set the diagonal to 1.0 */
  100.       Ai = A[i]; Ai[i] = 1.0;
  101.  
  102.       /* Divide by the pivot */
  103.       for ( j = 0; j < n; j++ )
  104.          {
  105.          Ai[j] *= Pivot;
  106.          }   /* for j */
  107.       Ai = B[i];
  108.       /* Divide by the pivot */
  109.       for ( j = 0; j < m; j++ )
  110.          {
  111.          Ai[j] *= Pivot;
  112.          }   /* for j */
  113.  
  114.       /* Loop through all the rows doing row
  115.       ** reduction */
  116.       for ( j = 0; j < n; j++ )
  117.          {
  118.          if ( j != i )
  119.             {
  120.             
  121.             /* Reuse the Pivot variable as a
  122.             ** dummy */
  123.             Ai = A[i]; Aj = A[j]; Pivot = Aj[i];
  124.             Aj[i] = 0.0;
  125.  
  126.             /* Row reduce the matrix */
  127.             for ( k = 0; k < n; k++ )
  128.                {
  129.                Aj[k] -= Ai[k] * Pivot;
  130.                }   /* for k */
  131.             Ai = B[i]; Aj = B[j];
  132.             /* Row reduce the solution vectors */
  133.             for ( k = 0; k < m; k++ )
  134.                {
  135.                Aj[k] -= Ai[k] * Pivot;
  136.                }   /* for k */
  137.  
  138.             }   /* if j */
  139.          }   /* for j */
  140.  
  141.       }   /* for i */
  142.  
  143.    return ( 0 );
  144.  
  145.    }   /* function gauss_jordan */
  146.  
  147. /* End of File */
  148.